Caustics and Intermittency in Turbulent Suspensions of Heavy Particles 
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The statistics of velocity differences between very heavy inertial particles suspended in an in- 
compressible turbulent flow is found to be extremely intermittent. When particles are separated 
by distances within the viscous subrange, the competition between quiet regular regions and multi- 
valued caustics leads to a quasi bi-fractal behavior of the particle velocity structure functions, with 
high-order moments bringing the statistical signature of caustics. Contrastingly, for particles sepa- 
rated by inertial-range distances, the velocity-difference statistics is characterized in terms of a local 
Holder exponent, which is a function of the scale-dependent particle Stokes number only. Results 
are supported by high-resolution direct numerical simulations. It is argued that these findings might 
have implications in the early stage of rain droplets formation in warm clouds. 

PACS numbers: 47.27.-i, 47.10.-g 



Two effects have recently been singled out to explain the 
speed-up of collisions between small finite-size particles 
suspended in turbulent flows [H, 0j S] ■ preferential con- 
centration, that is the development of strong inhomo- 
geneities in their spatial distribution (Fig. [Th.) 0, H, [|| , 
and the formation of fold caustics (also called the sling ef- 
fect), which results in large probabilities that close parti- 
cles have important velocity differences (Fig. [TJd) @, H, 0] • 
Improving the collision kernels used in kinetic models for 
atmospheric physics, astrophysics, and engineering re- 
quires quantifying precisely the individual contribution 
of these two effects and, in particular, to what extent 
turbulence might affect them, i.e. how they depend on 
the Taylor-scale Reynolds number R\ of the flow 1(| HI • 
In this Letter, to straighten out these questions, we 
consider suspensions of small, heavy, and dilute parti- 
cles, which is a setting relevant to the early stage of rain 
droplets formation in clouds Q. In these conditions, par- 
ticles are simply dragged by viscous forces and each in- 
dividual trajectory X(t) solves the equation 



between particles is evaluated using the ghost-particle 
approach [Bj], which assumes that particles can occupy 
any point of space independently of the positions of oth- 
ers. This approximation is valid in the asymptotics of 
very diluted suspensions, and has the advantage of rely- 
ing on stationary dynamical statistics: the geometrical 
collision rate is then determined by the value at r = 2a 
of the approaching rate 0] 



K,(r; St) 
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R = r and R < 0) p 2 (r) 



(2) 



Here R denotes the distance between two particles with 
Stokes number St, and p 2 its probability density. The 
average is performed over time and particle pairs, with 
the condition to be separated by a distance r and to ap- 
proach each other. Clearly the behavior of «(r; St) pre- 
scribes the dependence of the collision rate upon the par- 



tX = -X+u(X,t) 



(1) 



where dots denote time derivatives, u the fluid veloc- 
ity field, solution of the incompressible Navier-Stokes 
equation, and r = 2a 2 a/(9is) the Stokes time, depend- 
ing on particle radius a, density contrast a with the 




fluid, and fluid kinematic viscosity v, see [i^|J tor a re- 
cent review. The importance of inertia in the particle 
dynamics is quantified by the Stokes number St = t/t v , 
where r n is the fluid eddy turnover time associated to 
the Kolmogorov dissipative scale rj. The collision rate 



FIG. 1: (a) Snapshot of the position of particles for St = 2 
in a slice of size 5r/x 1007/ x IOO77 for R\ ps 400. (b) Particle 
velocity field in the same slice for a larger Stokes, St — 20, 
showing the existence of regions where particles have different 
velocities (highlighted by gray and black arrows). 
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tide attributes (size and mass density contrast). Caustics 
and preferential concentration (Fig. [I]) intricately appear 
in @ affecting the conditional average of the velocity 
difference R and the r-dependence of P2, respectively. In 
particular, as shown in 13], P2{r), which is straightfor- 
wardly related to the radial distribution function of [jj, 
behaves as a power law in the dissipative range, namely 
P2M oc r D2_1 for r<i), where D2 E [0 : 3] is the cor- 
relation dimension of the particle distribution and non- 
trivially depends on the Stokes number. 

We focus here on the velocity contribution by study- 
ing the behavior as a function of the separation r of the 
longitudinal particle velocity structure functions 



(3) 



S p {r;St) = (\R\P R = 



The choice of defining structure functions with abso- 
lute values is motivated by the definition of the colli- 
sion kernel (via the approaching rate), since we do not 
expect important differences between negative and pos- 
itive velocity fluctuations. One can therefore estimate: 
n(r) oc r^^S^St) (see Q). Evaluating S p (r;St) for 
values of p different from 1, besides providing a more 
complete characterization of the velocity statistics, al- 
lows one to account also for fluctuations of the local ap- 
proaching rate, which can vary significantly from place to 
place. In the limit of small inertia, the particle dynamics 
approaches that of tracers and consequently the velocity 
difference R is essentially coincident with the fluid longi- 
tudinal increment over a separation R. Conversely, when 
St — > co, particles move ballistically in the flow with un- 
corrected velocities and the structure functions S p (r; St) 
become independent of r. For intermediate values of the 
Stokes number, one expects a non-trivial behavior of S p 
as a function of r and St. Data analyzed in this study are 
from a direct numerical simulation at R\ « 400 described 
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Figure [2] represents the behavior of the second-order 
structure function 5*2 (r; Si), measured in direct numer- 
ical simulations, for two different values of the carrier 




FIG. 2: (color online) Second-order longitudinal velocity 
structure function for particles with various Stokes numbers 
and for two Reynolds numbers. 
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FIG. 3: (color online) Scaling exponent, £1, vs. St for two 
values of R\. Inset: correlation dimension D2 vs. St. 



flow Reynolds number (see [14| for details on the simu- 
lations). One distinguishes different regimes, depending 
whether r is within the dissipative or inertial range of 
the turbulent carrier flow. While the dissipative-range 
behavior directly relates to inter-particle collisions, the 
inertial-range behavior has important implications on the 
relative dispersion of particles in turbulent flow [3] in 
general and for pollution models in particular. In the 
sequel we investigate these two regimes separately. 

In the dissipative range, the structure functions display 
a power-law behavior S p (r; St) oc r^ p . The two asymp- 
totics of weak and strong inertia imply that £ p ~ p for 
St <C 1 and £ p — > for St — > 00. For intermediate values 
of the Stokes number, p 1— > £ p (St) is a convex function 
of the order p with values in [0 : p] . Figure [3] shows the 
first-order exponent £1 as a function of the Stokes num- 
ber. One can clearly observe that for St = 0(1), the 
exponent £1 takes non-trivial values spanning the whole 
interval [0:1]. The present accuracy of data does not al- 
low for settling either the issue of a possible saturation of 
the exponent to the limiting values at the two extrema, 
nor a possible dependence of the exponent upon R\. De- 
spite a factor two in R\, data differ by less than the 
errors made in the determination of the exponents or in 
the value of t„ that enters the definition of St. 

At first glance the continuous variation of the exponent 
£1 from 1 to at increasing St seems inconsistent with 
a naive picture of the role of caustics in velocity statis- 
tics. Fold caustics are a part of catastrophe theory [r|; 
they occur when fast particles catch up with slower ones 
to create regions where several velocities can be found at 
the same location as in Fig.QJ). If particles conserve their 
velocity and move ballistically, such caustics will extend 
over the whole domain (whence the analogy with caustics 
formed by light rays [8]). The typical velocity difference 
between two particles becomes in that case independent 
of their distance, meaning that structure functions tend 
to a constant as r — ► 0, and thus £ p = 0. However, 
there are two clear reasons why this continuous-field pic- 
ture may fail. First, because of their dissipative dynam- 
ics, particles concentrate on dynamical attractors in the 
position- velocity phase space Q- Such sets are fractal 
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FIG. 4: (color online) Scaling exponents £ p of the particle 
velocity structure functions S p for various St and R\ ~ 400. 
Inset: saturation exponents Coo as a function of St for two 
values of R\ . Exponents are obtained by measuring the mean 
logarithmic derivative of S p (r) in 0.2 < (r/rj) <2; errors corre- 
spond to the largest deviations observed in the fitting range. 

and correlated with the fluid and lead to the formation 
of caustics of various strength with non-trivial probabil- 
ities. Second, as the particle velocity relaxes to the fluid 
velocity, the spatio-temporal extent of such caustics may 
also have complex statistical properties. 

To better quantify the contribution from caustics, we 
extend our investigation to the particle velocity scal- 
ing exponents Cp's with orders p other than 1, shown 
in Fig. [4] for various values of the Stokes number. At 
small orders, the exponents are almost tangent to the 
line Cp = p while, at large orders, they approach or satu- 
rate to an asymptotic value Coo , monotonically decreasing 
with St as shown in the inset. Numerical data suggests 
that Coo oc \n(l/St) for St < 7 and that Coo - for 
St>7. The current numerical accuracy does not enable 
distinguishing between a sharp and smooth transition at 
St ~ 7. However, we notice that the estimated value 
of such a critical Stokes number is very close to that 
for which D 2 w 3 (see inset of Fig. [3]). As discussed in 
flpj ] a saturation of D2 to the space dimension is indeed 
expected for St values at which caustics become domi- 
nating. In this respect, we also notice that the satura- 
tion exponent Coo can be interpreted as the co-dimension 
of large fold caustics associated to order-unity velocity 
jumps. Indeed, such caustics contribute to the structure 
function S p (r; St) a term of the form V p P{r) where V p 
is the p-th order moment of the velocity difference in- 
side the caustics and P(r) is the probability of having 
such a caustic present in a box of size r. The satura- 
tion of the scaling exponents suggests that P(r) oc r^°° , 
so that D( c ) = 3 — Coo is the (statistical) Hausdorff di- 
mension of the set of caustics. At smaller orders, the 
statistics is dominated by other events for which one can 
figure out two conceivable scenarios. A first possibility 
is that caustics distribution spans all possible sizes with 
non-trivial co-dimensions, i.e. is a multi-fractal. In this 
case they affect all orders and give rise to multiscaling 
and to a non-trivial behavior of the exponents Cp before 
the saturation 17J. The alternative possibility is that 




FIG. 5: (color online) Probability density of the (rescaled) 
longitudinal velocity difference a = R/R for various values of 
r and for St = 3.3 and R\ ~ 185. Inset: same for the right 
tail in log-log coordinates. 

the caustics are randomly distributed with a typical size 
and dominate the velocity statistics at large moments 
only, while small orders are controlled by the smooth re- 
gions of the particle velocity. In that case the structure 
function would display a bi-fractal behavior similar to 
that present in random solutions to the Burgers equa- 
tion (see, e.g., [lij]), namely Cp = V for V < Coo and 
Cp = Coo for P > Coo- Current numerical results do not 
permit to distinguish between these two possibilities. As 
seen from Fig. |4j the measured exponents show some de- 
viations from the bi-fractal behavior. However as already 
observed in other settings , this apparent multiscaling 
could be an artifact due to the presence of sub-leading 
terms or logarithmic corrections. 

To further disentangle the question of the contribution 
of caustics to velocity scaling, we investigate the statis- 
tical properties of a = R/R, which can interpreted as 
a longitudinal velocity gradient of an effective particle 
velocity field. This quantity is at the center of much 
work devoted to the relative motion of a pair of parti- 
cles in time-uncorrelated flows [1, [5(1 21, There, 



the dynamics of a becomes independent of R at very 
small scales, a far from obvious feature for particles trans- 
ported by real flows, where time correlations and struc- 
tures play important roles. Further, results in random 
flows suggest that the conditional probability density 
p(a I R = r) is independent of a at small scales and has 
power-law tails. As seen from Fig. [5j numerical mea- 
surements in turbulent flows suggest features similar to 
those of structure-less random flows. The core of the 
distributions associated to different scales r collapse for 
H iS cr*{r) on a distribution with a fat, almost algebraic 
behavior. Interestingly, the associated power-law expo- 
nent is close to — (1 + Coo), suggesting that (<j p ) diverges 
for p > Coo, a behavior favoring the bi-fractal scenario. 



Indeed S p (r; St) = r P(a p \R = r} ~ 
p such that (a v ) < 00. However 



rP(a p ) for r -> and 
for |er| > a*(r), the 
distributions display stretched-exponential tails, whose 
contribution to the structure function is for the moment 
unsettled. They are connected to the caustics size prob- 
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FIG. 6: (color online) First-order local exponent £i(r) as a 
function of the local Stokes number St(r) for R\ ~ 400. The 
horizontal dashed line represents tracers K41 expectation. 



ability distribution and could lead to multiscaling. A 
related open question is the non-trivial entanglement be- 
tween clusters and large velocity differences, as already 
stressed in random flows (23[. This latter feature might 
imply an intricate dependence on r of velocity difference 
statistics, that might lead to multiscaling. Settling nu- 
merically the issue of bi- versus multiscaling would re- 
quire to explore systems were statistics can be handled 
in a more systematic way, as for instance in random cor- 
related flows or real flows at smaller Reynolds numbers. 

We finally turn to the behavior of the velocity struc- 
ture function for separations within the inertial range of 
turbulence, i.e. for r\ <C r <C L. As seen from Fig. [2j par- 
ticle velocity structure functions recover the fluid ones 
when r becomes very large. Indeed as r increases the 
associated eddy turnover time grows as r 2 / 3 (where we 
used the Kolmogorov 1941 scaling) so that the effective 
strength of inertia reduces. Similarly to random self- 
similar carrier flows [23], this effect can be put on a 
quantitative ground in terms of a scale-dependent Stokes 
number St(r) — e 1 / 3 T/r 2 / 3 defined as the ratio between 
the particle response time and the turnover time associ- 
ated to the scale r, where e denotes the mean dissipa- 
tion rate of kinetic energy. We check whether the local 
scaling exponent ^ P (r;r) = (dlnS , p (r; St))/ (din r) does 
depend on St(r) only, as observed in random self-similar 
flows [13] • Figure [5] shows a good collapse of the values of 
£i(r, t) associated to various r and of r, once represented 
as a function of St(r). Moreover, the curve £i(St(r)) has 
a shape qualitatively very similar to that of £,i(St) ob- 
served in the dissipative range and shown in Fig. [3l this 
fact is relevant to heavy particle dispersion in turbulent 
flows Let us stress that data corresponding to small 
St(r) in Fig.[6]show deviations from the K41 scaling that 
are similar to those expected for tracers-like statistcs. 

To conclude we briefly discuss the applicability of the 
present results to atmospheric physics. The main short- 
coming of the proposed approach is that the gravity force 
is neglected. As observed in 24J for dynamics of water 
droplets in warm clouds, gravitational settling is found 
to dominate the statistics of velocity differences between 



particles. However, this effect acts mainly between parti- 
cles of different sizes that fall at different speeds. Present 
results should apply to earlier stage of rain formation 
during which the droplets are almost mono-disperse and 
might play an important role in explaining the observed 
fast spectral broadening observed in clouds. 
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